Exploring Spectrum Sensitivies to Species Omission

For the sensitivity testing we estimated the size spectrum slope and intercept for every year and within each survey area an additional time, at each pass one of the species was omitted from the analysis. This was done to capture the relative importance each species had on the estimation of the biomass size spectra.

Loading Bio Data

The biological data used as a starting point is the same “node” used for the size spectrum estimations using the full community. This is loaded so that we can interrogate the overall biomass patterns that correspond to the results we see from the sensitivity testing.

# Access Biological Data with {gmRi}
# nmfs_bio <- gmRi::gmri_survdat_prep(survdat_source = "bio")
withr::with_dir(rprojroot::find_root('_targets.R'), 
                tar_load(nefsc_1g_binned))


# Get Yearly and Regionally averaged Biomasses
yr_szn_overall <- nefsc_1g_binned %>% 
  group_by(Year, survey_area) %>% 
  summarise(overall_bio = sum(strat_total_lwbio_s, na.rm = T),
            overall_abund = sum(strat_total_abund_s, na.rm = T),
            .groups = "drop")

# Get what it is for each species
sp_yr_szn <- nefsc_1g_binned %>% 
  group_by(comname, spec_class, fishery, hare_group, Year, survey_area) %>% 
  summarise(sp_bio = sum(strat_total_lwbio_s, na.rm = T),
            sp_abund = sum(strat_total_abund_s, na.rm = T),
            .groups = "drop")


# Join in the overall, get the fraction and what the percent would be
perc_bio <- left_join(sp_yr_szn, yr_szn_overall, by = c("Year", "survey_area")) %>% 
  mutate(
    bio_frac = sp_bio / overall_bio,
    abund_frac = sp_abund / overall_abund,
    bio_perc = bio_frac * 100,
    abund_perc = abund_frac * 100
  )

Overall Biomass Distributions

# perc_bio overall rankings
region_overalls <- perc_bio %>% 
  group_by(comname, spec_class, fishery, hare_group, survey_area) %>% 
  summarise(avg_bio_perc   = mean(bio_perc, na.rm = T),
            avg_abund_perc = mean(abund_perc, na.rm = T),
            .groups = "drop")

Biomass

# Treemaps
ggplot(region_overalls) +
  geom_treemap(aes(area = avg_bio_perc, fill = comname)) +
  geom_treemap_text(aes(area = avg_bio_perc, fill = comname, label = comname)) +
  facet_wrap(~survey_area, ncol = 2) +
  scale_fill_gmri() + 
  theme(legend.position = "none") +
  labs(title = "Fraction of Stratified Biomass: 1970-2019")

Abundance

# Treemaps
ggplot(region_overalls) +
  geom_treemap(aes(area = avg_abund_perc, fill = comname)) +
  geom_treemap_text(aes(area = avg_abund_perc, fill = comname, label = comname)) +
  facet_wrap(~survey_area, ncol = 2) +
  scale_fill_gmri() + 
  theme(legend.position = "none") +
  labs(title = "Fraction of Stratified Abundance: 1970-2019")

Yearly Patterns

Biomass

# Play around with plotting yearly influences
ggplot(perc_bio, aes(Year, bio_perc)) +
  # geom_point(aes(color = comname), show.legend = F) +
  # geom_line(aes(color = comname, group = comname), show.legend = F) +
  geom_col(aes(fill = spec_class)) +
  facet_grid(survey_area~.) +
  scale_color_gmri() +
  scale_fill_gmri() +
  theme(legend.position = "bottom") +
  labs(y = "Percent of Yearly Biomass",
       x = NULL,
       fill = "")

# 

Abundance

# Play around with plotting yearly influences
ggplot(perc_bio, aes(Year, abund_perc)) +
  # geom_point(aes(color = comname), show.legend = F) +
  # geom_line(aes(color = comname, group = comname), show.legend = F) +
  geom_col(aes(fill = spec_class)) +
  facet_grid(survey_area~.) +
  #facet_grid(spec_class~survey_area) +
  scale_color_gmri() +
  scale_fill_gmri() +
  theme(legend.position = "bottom") +
  labs(y = "Percent of Yearly Biomass",
       x = NULL,
       fill = "")

# 

Loading Sensitivity Results

# Access Biological Data with {gmRi}
# nmfs_bio <- gmRi::gmri_survdat_prep(survdat_source = "bio")
withr::with_dir(rprojroot::find_root('_targets.R'), 
                tar_load(species_sensitivity_shifts))

# Make year a numeric value
sens_results <- species_sensitivity_shifts %>%
  mutate(Year = as.numeric(Year)) 


# Join in the biomass info
sens_results <- sens_results %>% 
  select(-spec_class) %>% 
  left_join(perc_bio,
            by = c("Year" = "Year", 
                   "survey_area" = "survey_area", 
                   "spec_omit" = "comname"))

Single Species Removal Impacts

sens_ranks <- sens_results %>% 
  group_by(
    `Species Removed` = spec_omit, 
    Area = survey_area, 
    Guild = spec_class, 
    Fishery = fishery) %>% 
  summarise(`Avg. Slope Shift` = mean(slope_shift, na.rm = T),
            `Avg. Slope Shift (z)` = mean(slope_shift_z, na.rm = T),
            .groups = "drop") %>% 
  mutate(across(where(is.numeric), round, 2))

Slope Increase

# Slope Steepeners
sens_ranks %>% 
  arrange(`Avg. Slope Shift (z)`) %>% 
  head(8) %>% 
  gt::gt()  %>%
  gt::tab_header(
    title = "Top 8 Steepeners",
    subtitle = "Removing these species causes slope to steepen"
  )
Top 8 Steepeners
Removing these species causes slope to steepen
Species Removed Area Guild Fishery Avg. Slope Shift Avg. Slope Shift (z)
bullnose ray SNE NA Not Labelled -0.41 -2.16
silver hake GoM Groundfish Commercially Targeted -0.20 -1.08
silver hake SNE Groundfish Commercially Targeted -0.14 -0.75
barndoor skate GB Elasmobranch Commercially Targeted -0.11 -0.59
butterfish MAB Pelagic Not Commercially Targeted -0.08 -0.58
silver hake GB Groundfish Commercially Targeted -0.11 -0.56
roughtail stingray SNE NA Not Labelled -0.09 -0.49
spotted hake MAB Groundfish Not Commercially Targeted -0.06 -0.44

Slope Flattener

# Slope Flatteners
sens_ranks %>% 
  arrange(desc(`Avg. Slope Shift (z)`)) %>% 
  head(8) %>% 
  gt::gt()  %>%
  gt::tab_header(
    title = "Top 8 Flatteners",
    subtitle = "Removing these species causes slope to flatten"
  )
Top 8 Flatteners
Removing these species causes slope to flatten
Species Removed Area Guild Fishery Avg. Slope Shift Avg. Slope Shift (z)
spiny dogfish MAB Elasmobranch Commercially Targeted 0.13 0.92
spiny dogfish SNE Elasmobranch Commercially Targeted 0.15 0.80
winter skate GB Elasmobranch Commercially Targeted 0.04 0.22
spiny dogfish GoM Elasmobranch Commercially Targeted 0.04 0.21
atlantic cod GB Groundfish Commercially Targeted 0.04 0.19
spiny dogfish GB Elasmobranch Commercially Targeted 0.03 0.17
atlantic cod GoM Groundfish Commercially Targeted 0.02 0.12
sand tiger MAB NA Not Labelled 0.02 0.12

Regional Influence Heatmaps

# Do we need to show the relative influence?

# Things seem to wash out with the raw z-scores
# Maybe they should be adjusted/weighted for how much any single omission
# pulls a given year compared to either the 
# average z-shift for that year+region
# or by the max that it shifts....
# plot heatmap
influence_heatmap <- function(sens_results, 
                              filt_area, 
                              lims = c(-1, 1),
                              rescale_z = F){
  
  # Make breaks to indicate limits are truncating
  blabs <- seq(from = lims[1], to = lims[2], length.out = 5)
  blabs[c(1,5)] <- c(str_c(">", blabs[1]), 
                     str_c("<", blabs[5]))
  
  # Filter data
  filt_dat <- sens_results %>%
    mutate(Year = as.numeric(Year)) %>%
    filter(survey_area == filt_area)
  
  # Rescale?
  if(rescale_z == T){
    message("Rescaling not finished")
    
    
  }
                     
  
  # Make heatmap
  filt_dat %>%
  ggplot(aes(Year, fct_rev(spec_omit), fill = slope_shift_z)) +
    geom_tile() +
    scale_fill_distiller(palette = "RdBu",
                         limits = lims, 
                         labels = blabs,
                         oob = scales::oob_squish, 
                         na.value = "gray20") +
    scale_x_continuous(expand = expansion(add = c(0,0))) +
    facet_wrap(~survey_area, ncol = 2, scales = "free") +
    #facet_grid(spec_class~survey_area, scales = "free") +
    #facet_grid(fishery~survey_area, scales = "free") +
    theme(legend.position = "bottom", 
          legend.title = element_text(angle = 0), 
          axis.text.y = element_text(size = 7),
          panel.background = element_rect(fill = "gray40"),
          panel.grid = element_line(color = "transparent")) +
    guides(fill = guide_colorbar(title.position = "top", 
                                 title.hjust = 0.5, 
                                 barwidth = unit(3.25, "in"))) +
    labs(y = "Species omitted", 
         fill = "Shift in Slope (z-score)")
  
  
}

Gulf of Maine

# GOM
sens_results %>% 
  #filter(spec_class == "Groundfish") %>% 
influence_heatmap(filt_area =  "GoM", lims = c(-2,2))

Georges Bank

# GOM
influence_heatmap(sens_results, "GB", lims = c(-2,2))

Southern New England

# GOM
influence_heatmap(sens_results, "SNE", lims = c(-4,4))

Mid-Atlantic Bight

# GOM
influence_heatmap(sens_results, "MAB", lims = c(-4,4))

Correlation to Biomass Fraction

sens_results %>% 
  ggplot(aes(bio_frac, slope_shift_z)) +
  geom_point(alpha = 0.3) +
  geom_smooth(formula = "y ~ x", method = "lm") +
  labs(x = "Fration of Biomass", y = "Shift in Spectra Slope (z)")

Regional Influence Horizon Plots

I think this might help get around the “wash out” of the heatmap. The color scale should be longer, or potentially just not center on white…

Singe Species Checks:

# plot a species:
# multiple panel plot function
# Function for plotting a species that stands out
plot_spectra_influence <- function(sens_results, filter_spec, filter_area){

  # Grab the relevant data
  filt_d <- sens_results %>%
    filter(spec_omit == filter_spec,
           survey_area == filter_area)
  
  # Make Plot Label
  p_label <- str_c(str_to_title(filter_spec), " Omission")
  
  
  # 1. Actual Value Shift
  # 1a. Slope
  splot <- filt_d %>%
    ggplot(aes(Year, l10_slope_strat)) +
    geom_line(aes(color = "All Species"), group = 1) +
    geom_point(aes(color = "All Species")) +
    geom_line(aes(y = omit_slope_strat, color = p_label), group = 1) +
    geom_point(aes(y = omit_slope_strat, color = p_label)) +
    facet_wrap(~survey_area, ncol = 1) +
    scale_color_gmri() +
    theme(legend.position = "bottom") +
    labs(color = NULL,
         y = "Biomass Spectrum Slope",
         x = NULL)
  
  # 2. What the z-score shift is
  z_shift <- filt_d %>% 
    ggplot(aes(Year, slope_shift_z)) +
    geom_line() +
    geom_point() +
    labs(y = "Slope Shift (z)",
         x = NULL)
  
  # Biomass Fraction Each Year:
  bfrac <- filt_d %>%
    ggplot(aes(Year, bio_frac)) +
    geom_line() +
    geom_point() +
    scale_y_continuous(labels = scales::percent) +
    labs(y = "Biomass Fraction")
  
  
  # Do they fit
  plotsemble <- splot / z_shift / bfrac
  return(plotsemble)
}

Haddock (nowhere)

plot_spectra_influence(sens_results = sens_results, 
                       filter_spec = "haddock", 
                       filter_area = "GoM")

Spiny Dogfish

plot_spectra_influence(sens_results = sens_results, 
                       filter_spec = "atlantic herring", 
                       filter_area = "SNE")

Silver Hake (GoM & GB)

plot_spectra_influence(sens_results = sens_results, 
                       filter_spec = "silver hake", 
                       filter_area = "GoM")

plot_spectra_influence(sens_results = sens_results, 
                       filter_spec = "silver hake", 
                       filter_area = "GB")

Butterfish

plot_spectra_influence(sens_results = sens_results, 
                       filter_spec = "butterfish", 
                       filter_area = "GoM")

Atlantic Cod

plot_spectra_influence(sens_results = sens_results, 
                       filter_spec = "atlantic cod", 
                       filter_area = "GoM")